MATH50003 Numerical Analysis Lecture Notes

Author: Dr. Sheehan Olver | Github Source

Asymptotics and Computational Cost

YouTube Lectures: Asymptotics and Computational Cost

We introduce Big-O, little-o and asymptotic notation and see how they can be used to describe computational cost.

1. Asymptotics as $n\to \infty$$n → ∞$
2. Asymptotics as $x\to x_0$$x → x_0$
3. Computational cost

1. Asymptotics as $n\to \infty$$n → ∞$

Big-O, little-o, and "asymptotic to" are used to describe behaviour of functions at infinity.

Definition (Big-O)

$$f(n)=O(\varphi (n))\text{(as }n\to \infty \text{)}$$

means

$$\left|\frac{f(n)}{\varphi (n)}\right|$$

is bounded for sufficiently large $n$$n$. That is, there exist constants $C$$C$ and $N_0$$N_0$ such that, for all $n\ge N_0$$n \geq N_0$, $|\frac{f(n)}{\varphi (n)}|\le C$$|{f(n) \over ϕ(n)}| \leq C$.

Definition (little-O)

$$f(n)=o(\varphi (n))\text{(as }n\to \infty \text{)}$$

means

$$\underset{n\to \infty }{lim}\frac{f(n)}{\varphi (n)}=0.$$

Definition (asymptotic to)

$$f(n)\sim \varphi (n)\text{(as }n\to \infty \text{)}$$

means

$$\underset{n\to \infty }{lim}\frac{f(n)}{\varphi (n)}=1.$$

Example (integer in binary)

as

for $n\ge N_0=2$$n \geq N_0 = 2$.

as

as

Note we sometimes write $f(O(\varphi (n)))$$f(O(ϕ(n)))$ for a function of the form $f(g(n))$$f(g(n))$ such that $g(n)=O(\varphi (n))$$g(n) = O(ϕ(n))$.

Rules

We have some simple algebraic rules:

Proposition (Big-O rules)

$$\begin{array}{r}O(\varphi (n))O(\psi (n))=O(\varphi (n)\psi (n))\text{(as }n\to \infty \text{)}\\ O(\varphi (n))+O(\psi (n))=O(|\varphi (n)|+|\psi (n)|)\text{(as }n\to \infty \text{)}.\end{array}$$

2. Asymptotics as $x\to x_0$$x → x_0$

We also have Big-O, little-o and "asymptotic to" at a point:

Definition (Big-O)

$$f(x)=O(\varphi (x))\text{(as }x\to x_0\text{)}$$

means

$$\frac{|f(x)}{\varphi (x)|}$$

is bounded in a neighbourhood of $x_0$$x_0$. That is, there exist constants $C$$C$ and $r$$r$ such that, for all $0\le |x-x_0|\le r$$0 \leq |x - x_0| \leq r$, $|\frac{f(x)}{\varphi (x)}|\le C$$|{f(x) \over ϕ(x)}| \leq C$.

Definition (little-O)

$$f(x)=o(\varphi (x))\text{(as }x\to x_0\text{)}$$

means

$$\underset{x\to x_0}{lim}\frac{f(x)}{\varphi (x)}=0.$$

Definition (asymptotic to)

$$f(x)\sim \varphi (x)\text{(as }x\to x_0\text{)}$$

means

$$\underset{x\to x_0}{lim}\frac{f(x)}{\varphi (x)}=1.$$

Example (integer in binary)

Since

for some $t\in [0,x]$$t \in [0,x]$ and

provided $x\le 1$$x \leq 1$.

3. Computational cost

We will use Big-O notation to describe the computational cost of algorithms. Consider the following simple sum

which we might implement as:

xxxxxxxxxx
function sumsq(x)
    n = length(x)
    ret = 0.0
    for k = 1:n
        ret = ret + x[k]^2
    end
    ret
end

n = 100
x = randn(n)
sumsq(x)

Each step of this algorithm consists of one memory look-up (z = x[k]), one multiplication (w = z*z) and one addition (ret = ret + w). We will ignore the memory look-up in the following discussion. The number of CPU operations per step is therefore 2 (the addition and multiplication). Thus the total number of CPU operations is $2n$$2n$. But the constant $2$$2$ here is misleading: we didn't count the memory look-up, thus it is more sensible to just talk about the asymptotic complexity, that is, the computational cost is $O(n)$$O(n)$.

Now consider a double sum like:

which we might implement as:

xxxxxxxxxx
function sumsq2(x)
    n = length(x)
    ret = 0.0
    for k = 1:n
        for j = 1:k
            ret = ret + x[j]^2
        end
    end
    ret
end

n = 100
x = randn(n)
sumsq2(x)

Now the inner loop is $O(1)$$O(1)$ operations (we don't try to count the precise number), which we do $k$$k$ times for $O(k)$$O(k)$ operations as $k\to \infty$$k → ∞$. The outer loop therefore takes

operations.

Numbers

YouTube Lectures: Integers Floating Point Numbers Rounding Bouding Rounding Errors

Reference: Overton

In this chapter, we introduce the Two's-complement storage for integers and the IEEE Standard for Floating-Point Arithmetic. There are many possible ways of representing real numbers on a computer, as well as the precise behaviour of operations such as addition, multiplication, etc. Before the 1980s each processor had potentially a different representation for real numbers, as well as different behaviour for operations.
IEEE introduced in 1985 was a means to standardise this across processors so that algorithms would produce consistent and reliable results.

This chapter may seem very low level for a mathematics course but there are two important reasons to understand the behaviour of integers and floating-point numbers:

1. Integer arithmetic can suddenly start giving wrong negative answers when numbers become large.
2. Floating-point arithmetic is very precisely defined, and can even be used in rigorous computations as we shall see in the problem sheets. But it is not exact and its important to understand how errors in computations can accumulate.
3. Failure to understand floating-point arithmetic can cause catastrophic issues in practice, with the extreme example being the explosion of the Ariane 5 rocket.

In this chapter we discuss the following:

1. Binary representation: Any real number can be represented in binary, that is, by an infinite sequence of 0s and 1s (bits). We review binary representation.
2. Integers: There are multiple ways of representing integers on a computer. We discuss the the different types of integers and their representation as bits, and how arithmetic operations behave like modular arithmetic. As an advanced topic we discuss BigInt, which uses variable bit length storage.
3. Floating-point numbers: Real numbers are stored on a computer with a finite number of bits. There are three types of floating-point numbers: normal numbers, subnormal numbers, and special numbers.
4. Arithmetic: Arithmetic operations in floating-point are exact up to rounding, and how the rounding mode can be set. This allows us to bound errors computations.
5. High-precision floating-point numbers: As an advanced topic, we discuss how the precision of floating-point arithmetic can be increased arbitrary using BigFloat.

Before we begin, we load two external packages. SetRounding.jl allows us to set the rounding mode of floating-point arithmetic. ColorBitstring.jl implements functions printbits (and printlnbits) which print the bits (and with a newline) of floating-point numbers in colour.

xxxxxxxxxx
using SetRounding, ColorBitstring

1. Binary representation

Any integer can be presented in binary format, that is, a sequence of 0s and 1s.

Definition For $B_0,\dots ,B_p\in \{0,1\}$$B_0,\ldots,B_p \in \{0,1\}$ denote a non-negative integer in binary format by:

$$(B_p\dots B_1{B}_0{)}_2:=2^p{B}_p+\cdots +2B_1+B_0$$

For $b_1,b_2,\dots \in \{0,1\}$$b_1,b_2,\ldots \in \{0,1\}$, Denote a non-negative real number in binary format by:

$$(B_p\dots B_0.b_1{b}_2{b}_3\dots {)}_2=(B_p\dots B_0{)}_2+\frac{b_1}{2}+\frac{b_2}{2^2}+\frac{b_3}{2^3}+\cdots$$

First we show some examples of verifying a numbers binary representation:

Example (integer in binary) A simple integer example is $5=2^2+2^0=(101{)}_2$$5 = 2^2 + 2^0 = (101)_2$.

Example (rational in binary) Consider the number 1/3. In decimal recall that:

We will see that in binary

Both results can be proven using the geometric series:

provided $|z|<1$$|z| < 1$. That is, with $z=\frac{1}{4}$$z = {1 \over 4}$ we verify the binary expansion:

A similar argument with $z=1/10$$z = 1/10$ shows the decimal case.

2. Integers

On a computer one typically represents integers by a finite number of $p$$p$ bits, with $2^p$$2^p$ possible combinations of 0s and 1s. For unsigned integers (non-negative integers) these bits are just the first $p$$p$ binary digits: $(B_{p-1}\dots B_1{B}_0{)}_2$$(B_{p-1}\ldots B_1B_0)_2$.

Integers on a computer follow modular arithmetic:

Definition (ring of integers modulo $m$$m$) Denote the ring

$${\mathbb{Z}}_m:=\{0(\mathrm{m}\mathrm{o}\mathrm{d}m),1(\mathrm{m}\mathrm{o}\mathrm{d}m),\dots ,m-1(\mathrm{m}\mathrm{o}\mathrm{d}m)\}$$

Integers represented with $p$$p$-bits on a computer actually represent elements of ${\mathbb{Z}}_{2^p}$${\mathbb Z}_{2^p}$ and integer arithmetic on a computer is equivalent to arithmetic modulo $2^p$$2^p$.

Example (addition of 8-bit unsigned integers) Consider the addition of two 8-bit numbers:

The result is impossible to store in just 8-bits! It is way too slow for a computer to increase the number of bits, or to throw an error (checks are slow). So instead it treats the integers as elements of ${\mathbb{Z}}_{256}$${\mathbb Z}_{256}$:

We can see this in Julia:

xxxxxxxxxx
x = UInt8(255)
y = UInt8(1)
printbits(x); println(" + "); printbits(y); println(" = ")
printbits(x + y)

Example (multiplication of 8-bit unsigned integers) Multiplication works similarly: for example,

We can see this behaviour in code by printing the bits:

xxxxxxxxxx
x = UInt8(254) # 254 represented in 8-bits as an unsigned integer
y = UInt8(2)   #   2 represented in 8-bits as an unsigned integer
printbits(x); println(" * "); printbits(y); println(" = ")
printbits(x * y)

Signed integer

Signed integers use the Two's complemement convention. The convention is if the first bit is 1 then the number is negative: the number $2^p-y$$2^p - y$ is interpreted as $-y$$-y$. Thus for $p=8$$p = 8$ we are interpreting $2^7$$2^7$ through $2^8-1$$2^8-1$ as negative numbers.

Example (converting bits to signed integers) What 8-bit integer has the bits 01001001? Adding the corresponding decimal places we get:

xxxxxxxxxx
2^0 + 2^3 + 2^6

What 8-bit (signed) integer has the bits 11001001? Because the first bit is 1 we know it's a negative number, hence we need to sum the bits but then subtract 2^p:

xxxxxxxxxx
2^0 + 2^3 + 2^6 + 2^7 - 2^8

We can check the results using printbits:

xxxxxxxxxx
printlnbits(Int8(73))
printbits(-Int8(55))

Arithmetic works precisely the same for signed and unsigned integers.

Example (addition of 8-bit integers) Consider (-1) + 1 in 8-bit arithmetic. The number $-1$$-1$ has the same bits as $2^8-1=255$$2^8 - 1 = 255$. Thus this is equivalent to the previous question and we get the correct result of 0. In other words:

Example (multiplication of 8-bit integers) Consider (-2) * 2. $-2$$-2$ has the same bits as $2^{256}-2=254$$2^{256} - 2 = 254$ and $-4$$-4$ has the same bits as $2^{256}-4=252$$2^{256}-4 = 252$, and hence from the previous example we get the correct result of -4. In other words:

Example (overflow) We can find the largest and smallest instances of a type using typemax and typemin:

xxxxxxxxxx
printlnbits(typemax(Int8)) # 2^7-1 = 127
printbits(typemin(Int8)) # -2^7 = -128

As explained, due to modular arithmetic, when we add 1 to the largest 8-bit integer we get the smallest:

xxxxxxxxxx
typemax(Int8) + Int8(1) # returns typemin(Int8)

This behaviour is often not desired and is known as overflow, and one must be wary of using integers close to their largest value.

Variable bit representation (advanced)

An alternative representation for integers uses a variable number of bits, with the advantage of avoiding overflow but with the disadvantage of a substantial speed penalty. In Julia these are BigInts, which we can create by calling big on an integer:

xxxxxxxxxx
x = typemax(Int64) + big(1) # Too big to be an `Int64`

Note in this case addition automatically promotes an Int64 to a BigInt. We can create very large numbers using BigInt:

xxxxxxxxxx
x^100

Note the number of bits is not fixed, the larger the number, the more bits required to represent it, so while overflow is impossible, it is possible to run out of memory if a number is astronomically large: go ahead and try x^x (at your own risk).

Division

In addition to +, -, and * we have integer division ÷, which rounds down:

xxxxxxxxxx
5 ÷ 2 # equivalent to div(5,2)

Standard division / (or \ for division on the right) creates a floating-point number, which will be discussed shortly:

xxxxxxxxxx
5 / 2 # alternatively 2 \ 5

We can also create rational numbers using //:

xxxxxxxxxx
(1//2) + (3//4)

Rational arithmetic often leads to overflow so it is often best to combine big with rationals:

xxxxxxxxxx
big(102324)//132413023 + 23434545//4243061 + 23434545//42430534435

3. Floating-point numbers

Floating-point numbers are a subset of real numbers that are representable using a fixed number of bits.

Definition (floating-point numbers) Given integers $\sigma$$σ$ (the "exponential shift") $Q$$Q$ (the number of exponent bits) and $S$$S$ (the precision), define the set of Floating-point numbers by dividing into normal, sub-normal, and special number subsets:

$$F_{\sigma ,Q,S}:=F_{\sigma ,Q,S}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}}\cup F_{\sigma ,Q,S}^{\mathrm{s}\mathrm{u}\mathrm{b}}\cup F^{\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{a}\mathrm{l}}.$$

The normal numbers $F_{\sigma ,Q,S}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}}\subset \mathbb{R}$$F^{\rm normal}_{σ,Q,S} \subset {\mathbb R}$ are defined by

$$F_{\sigma ,Q,S}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}}=\{\pm 2^{q-\sigma }\times (1.b_1{b}_2{b}_3\dots b_S{)}_2:1\le q<2^Q-1\}.$$

The sub-normal numbers $F_{\sigma ,Q,S}^{\mathrm{s}\mathrm{u}\mathrm{b}}\subset \mathbb{R}$$F^{\rm sub}_{σ,Q,S} \subset {\mathbb R}$ are defined as

$$F_{\sigma ,Q,S}^{\mathrm{s}\mathrm{u}\mathrm{b}}=\{\pm 2^{1-\sigma }\times (0.b_1{b}_2{b}_3\dots b_S{)}_2\}.$$

The special numbers $F^{\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{a}\mathrm{l}}\not\subset \mathbb{R}$$F^{\rm special} \not\subset {\mathbb R}$ are defined later.

Note this set of real numbers has no nice algebraic structure: it is not closed under addition, subtraction, etc. We will therefore need to define approximate versions of algebraic operations later.

Floating-point numbers are stored in $1+Q+S$$1 + Q + S$ total number of bits, in the format

The first bit ($s$$s$) is the sign bit: 0 means positive and 1 means negative. The bits $q_{Q-1}\dots q_0$$q_{Q-1}\ldots q_0$ are the exponent bits: they are the binary digits of the unsigned integer $q$$q$:

Finally, the bits $b_1\dots b_S$$b_1\ldots b_S$ are the significand bits. If $1\le q<2^Q-1$$1 \leq q < 2^Q-1$ then the bits represent the normal number

If $q=0$$q = 0$ (i.e. all bits are 0) then the bits represent the sub-normal number

If $q=2^Q-1$$q = 2^Q-1$ (i.e. all bits are 1) then the bits represent a special number, discussed later.

IEEE floating-point numbers

Definition (IEEE floating-point numbers) IEEE has 3 standard floating-point formats: 16-bit (half precision), 32-bit (single precision) and 64-bit (double precision) defined by:

$$\begin{array}{rl}{F}_{16}& :=F_{15,5,10}\\ F_{32}& :=F_{127,8,23}\\ F_{64}& :=F_{1023,11,52}\end{array}$$

In Julia these correspond to 3 different floating-point types:

1. Float64 is a type representing double precision ($F_{64}$$F_{64}$). We can create a Float64 by including a decimal point when writing the number: 1.0 is a Float64. Float64 is the default format for scientific computing (on the Floating-Point Unit, FPU).
2. Float32 is a type representing single precision ($F_{32}$$F_{32}$). We can create a Float32 by including a f0 when writing the number: 1f0 is a Float32. Float32 is generally the default format for graphics (on the Graphics Processing Unit, GPU), as the difference between 32 bits and 64 bits is indistinguishable to the eye in visualisation, and more data can be fit into a GPU's limitted memory.
3. Float16 is a type representing half-precision ($F_{16}$$F_{16}$). It is important in machine learning where one wants to maximise the amount of data and high accuracy is not necessarily helpful.

Example (rational in Float32) How is the number $1/3$$1/3$ stored in Float32? Recall that

and since $125=(1111101{)}_2$$125 = (1111101)_2$ the exponent bits are 01111101. . For the significand we round the last bit to the nearest element of $F_{32}$$F_{32}$, (this is explained in detail in the section on rounding), so we have

and the significand bits are 01010101010101010101011. Thus the Float32 bits for $1/3$$1/3$ are:

xxxxxxxxxx
printbits(1f0/3)

For sub-normal numbers, the simplest example is zero, which has $q=0$$q=0$ and all significand bits zero:

xxxxxxxxxx
printbits(0.0)

Unlike integers, we also have a negative zero:

xxxxxxxxxx
printbits(-0.0)

This is treated as identical to 0.0 (except for degenerate operations as explained in special numbers).

Special normal numbers

When dealing with normal numbers there are some important constants that we will use to bound errors.

Definition (machine epsilon/smallest positive normal number/largest normal number) Machine epsilon is denoted

$${ϵ}_{\mathrm{m},S}:=2^{-S}.$$

When $S$$S$ is implied by context we use the notation ${ϵ}_{\mathrm{m}}$$ϵ_{\rm m}$. The smallest positive normal number is $q=1$$q = 1$ and $b_k$$b_k$ all zero:

$$min|F_{\sigma ,Q,S}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}}|=2^{1-\sigma }$$

where $|A|:=\{|x|:x\in A\}$$|A| := \{|x| : x \in A \}$. The largest (positive) normal number is

$$max{F}_{\sigma ,Q,S}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}}=2^{2^Q-2-\sigma }(1.11\dots 1{)}_2=2^{2^Q-2-\sigma }(2-{ϵ}_{\mathrm{m}})$$

We confirm the simple bit representations:

xxxxxxxxxx
σ,Q,S = 127,23,8 # Float32
εₘ = 2.0^(-S)
printlnbits(Float32(2.0^(1-σ))) # smallest positive Float32
printlnbits(Float32(2.0^(2^Q-2-σ) * (2-εₘ))) # largest Float32

For a given floating-point type, we can find these constants using the following functions:

xxxxxxxxxx
eps(Float32),floatmin(Float32),floatmax(Float32)

Example (creating a sub-normal number) If we divide the smallest normal number by two, we get a subnormal number:

xxxxxxxxxx
mn = floatmin(Float32) # smallest normal Float32
printlnbits(mn)
printbits(mn/2)

Can you explain the bits?

Special numbers

The special numbers extend the real line by adding $\pm \mathrm{\infty }$$\pm \infty$ but also a notion of "not-a-number".

Definition (not a number) Let $\mathrm{N}\mathrm{a}\mathrm{N}$${\rm NaN}$ represent "not a number" and define

$$F^{\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{a}\mathrm{l}}:=\{\mathrm{\infty },-\mathrm{\infty },\mathrm{N}\mathrm{a}\mathrm{N}\}$$

Whenever the bits of $q$$q$ of a floating-point number are all 1 then they represent an element of $F^{\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{a}\mathrm{l}}$$F^{\rm special}$. If all $b_k=0$$b_k=0$, then the number represents either $\pm \mathrm{\infty }$$\pm\infty$, called Inf and -Inf for 64-bit floating-point numbers (or Inf16, Inf32 for 16-bit and 32-bit, respectively):

xxxxxxxxxx
printlnbits(Inf16)
printbits(-Inf16)

All other special floating-point numbers represent $\mathrm{N}\mathrm{a}\mathrm{N}$${\rm NaN}$. One particular representation of $\mathrm{N}\mathrm{a}\mathrm{N}$${\rm NaN}$ is denoted by NaN for 64-bit floating-point numbers (or NaN16, NaN32 for 16-bit and 32-bit, respectively):

xxxxxxxxxx
printbits(NaN16)

These are needed for undefined algebraic operations such as:

xxxxxxxxxx
0/0

Example (many NaNs) What happens if we change some other $b_k$$b_k$ to be nonzero? We can create bits as a string and see:

xxxxxxxxxx
i = parse(UInt16, "0111110000010001"; base=2)
reinterpret(Float16, i)

Thus, there are more than one NaNs on a computer.

4. Arithmetic

Arithmetic operations on floating-point numbers are exact up to rounding. There are three basic rounding strategies: round up/down/nearest. Mathematically we introduce a function to capture the notion of rounding:

Definition (rounding) ${\mathrm{f}\mathrm{l}}_{\sigma ,Q,S}^{\mathrm{u}\mathrm{p}}:\mathbb{R}\to F_{\sigma ,Q,S}$${\rm fl}^{\rm up}_{σ,Q,S} : \mathbb R \rightarrow F_{σ,Q,S}$ denotes the function that rounds a real number up to the nearest floating-point number that is greater or equal. ${\mathrm{f}\mathrm{l}}_{\sigma ,Q,S}^{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}}:\mathbb{R}\to F_{\sigma ,Q,S}$${\rm fl}^{\rm down}_{σ,Q,S} : \mathbb R \rightarrow F_{σ,Q,S}$ denotes the function that rounds a real number down to the nearest floating-point number that is greater or equal. ${\mathrm{f}\mathrm{l}}_{\sigma ,Q,S}^{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}:\mathbb{R}\to F_{\sigma ,Q,S}$${\rm fl}^{\rm nearest}_{σ,Q,S} : \mathbb R \rightarrow F_{σ,Q,S}$ denotes the function that rounds a real number to the nearest floating-point number. In case of a tie, it returns the floating-point number whose least significant bit is equal to zero. We use the notation $\mathrm{f}\mathrm{l}$${\rm fl}$ when $\sigma ,Q,S$$σ,Q,S$ and the rounding mode are implied by context, with ${\mathrm{f}\mathrm{l}}^{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}$${\rm fl}^{\rm nearest}$ being the default rounding mode.

In Julia, the rounding mode is specified by tags RoundUp, RoundDown, and RoundNearest. (There are also more exotic rounding strategies RoundToZero, RoundNearestTiesAway and RoundNearestTiesUp that we won't use.)

WARNING (rounding performance, advanced) These rounding modes are part of the FPU instruction set so will be (roughly) equally fast as the default, RoundNearest. Unfortunately, changing the rounding mode is expensive, and is not thread-safe.

Let's try rounding a Float64 to a Float32.

xxxxxxxxxx
printlnbits(1/3)  # 64 bits
printbits(Float32(1/3))  # round to nearest 32-bit

The default rounding mode can be changed:

xxxxxxxxxx
printbits(Float32(1/3,RoundDown) )

Or alternatively we can change the rounding mode for a chunk of code using setrounding. The following computes upper and lower bounds for /:

xxxxxxxxxx
x = 1f0
setrounding(Float32, RoundDown) do
    x/3
end,
setrounding(Float32, RoundUp) do
    x/3
end

WARNING (compiled constants, advanced): Why did we first create a variable x instead of typing 1f0/3? This is due to a very subtle issue where the compiler is too clever for it's own good: it recognises 1f0/3 can be computed at compile time, but failed to recognise the rounding mode was changed.

In IEEE arithmetic, the arithmetic operations +, -, *, / are defined by the property that they are exact up to rounding. Mathematically we denote these operations as follows:

Note also that ^ and sqrt are similarly exact up to rounding.

Example (decimal is not exact) 1.1+0.1 gives a different result than 1.2:

xxxxxxxxxx
x = 1.1
y = 0.1
x + y - 1.2 # Not Zero?!?

This is because $\mathrm{f}\mathrm{l}(1.1)\ne 1+1/10$${\rm fl}(1.1) \neq 1+1/10$, but rather:

WARNING (non-associative) (non-associative) These operations are not associative! E.g. $(x\oplus y)\oplus z$$(x \oplus y) \oplus z$ is not necessarily equal to $x\oplus (y\oplus z)$$x \oplus (y \oplus z)$. Commutativity is preserved, at least. Here is a surprising example of non-associativity:

xxxxxxxxxx
(1.1 + 1.2) + 1.3, 1.1 + (1.2 + 1.3)

Can you explain this in terms of bits?

Bounding errors in floating point arithmetic

Before we dicuss bounds on errors, we need to talk about the two notions of errors:

Definition (absolute/relative error) If $\tilde{x}=x+{\delta }_{rma}=x(1+{\delta }_{\mathrm{r}})$$\tilde x = x + δ_{rm a} = x (1 + δ_{\rm r})$ then $|{\delta }_{\mathrm{a}}|$$|δ_{\rm a}|$ is called the absolute error and $|{\delta }_{\mathrm{r}}|$$|δ_{\rm r}|$ is called the relative error in approximating $x$$x$ by $\tilde{x}$$\tilde x$.

We can bound the error of basic arithmetic operations in terms of machine epsilon, provided a real number is close to a normal number:

Definition (normalised range) The normalised range ${\mathcal{N}}_{\sigma ,Q,S}\subset \mathbb{R}$${\cal N}_{σ,Q,S} \subset {\mathbb R}$ is the subset of real numbers that lies between the smallest and largest normal floating-point number:

$${\mathcal{N}}_{\sigma ,Q,S}:=\{x:min|F_{\sigma ,Q,S}|\le |x|\le max{F}_{\sigma ,Q,S}\}$$

When $\sigma ,Q,S$$σ,Q,S$ are implied by context we use the notation $\mathcal{N}$${\cal N}$.

We can use machine epsilon to determine bounds on rounding:

Proposition (rounding arithmetic) If $x\in \mathcal{N}$$x \in {\cal N}$ then

$${\mathrm{f}\mathrm{l}}^{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}}(x)=x(1+{\delta }_x^{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}})$$

where the relative error is

$$\begin{array}{rl}|{\delta }_x^{\mathrm{n}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}}|& \le \frac{{ϵ}_{\mathrm{m}}}{2}\\ |{\delta }_x^{\mathrm{u}\mathrm{p}/\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}}|& <{ϵ}_{\mathrm{m}}.\end{array}$$

This immediately implies relative error bounds on all IEEE arithmetic operations, e.g., if $x+y\in \mathcal{N}$$x+y \in {\cal N}$ then we have

where (assuming the default nearest rounding) $|{\delta }_1|\le \frac{{ϵ}_{\mathrm{m}}}{2}.$$|\delta_1| \leq {ϵ_{\rm m} \over 2}.$

Example (bounding a simple computation) We show how to bound the error in computing

using floating-point arithmetic. First note that 1.1 on a computer is in fact $\mathrm{f}\mathrm{l}(1.1)$${\rm fl}(1.1)$. Thus this computation becomes

First we find

where (note ${\delta }_1{\delta }_3$$δ_1 δ_3$ and ${\delta }_2{\delta }_3$$δ_2 δ_3$ are tiny so we just round up our bound to the nearest decimal)

Thus the computation becomes

where the absolute error is

Indeed, this bound is bigger than the observed error:

xxxxxxxxxx
abs(3.6 - (1.1+1.2+1.3)), 4.8eps()

Arithmetic and special numbers

Arithmetic works differently on Inf and NaN and for undefined operations. In particular we have:

xxxxxxxxxx
1/0.0        #  Inf
1/(-0.0)     # -Inf
0.0/0.0      #  NaN
  
Inf*0        #  NaN
Inf+5        #  Inf
(-1)*Inf     # -Inf
1/Inf        #  0.0
1/(-Inf)     # -0.0
Inf - Inf    #  NaN
Inf ==  Inf  #  true
Inf == -Inf  #  false

NaN*0        #  NaN
NaN+5        #  NaN
1/NaN        #  NaN
NaN == NaN   #  false
NaN != NaN   #  true

Special functions (advanced)

Other special functions like cos, sin, exp, etc. are not part of the IEEE standard. Instead, they are implemented by composing the basic arithmetic operations, which accumulate errors. Fortunately many are designed to have relative accuracy, that is, s = sin(x) (that is, the Julia implementation of $\sin x$$\sin x$) satisfies

where $|\delta |<c{ϵ}_{\mathrm{m}}$$|\delta| < cϵ_{\rm m}$ for a reasonably small $c>0$$c > 0$, provided that $x\in {\mathrm{F}}^{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}}$$x \in {\rm F}^{\rm normal}$. Note these special functions are written in (advanced) Julia code, for example, sin. WARNING (sin(fl(x)) is not always close to sin(x)) This is possibly a misleading statement when one thinks of $x$$x$ as a real number. Consider $x=\pi$$x = \pi$ so that $\sin x=0$$\sin x = 0$. However, as $\mathrm{f}\mathrm{l}(\pi )\ne \pi$${\rm fl}(\pi) \neq \pi$. Thus we only have relative accuracy compared to the floating point approximation:

xxxxxxxxxx
π₆₄ = Float64(π)
πᵦ = big(π₆₄) # Convert 64-bit approximation of π to higher precision. Note its the same number.
abs(sin(π₆₄)), abs(sin(π₆₄) - sin(πᵦ)) # only has relative accuracy compared to sin(πᵦ), not sin(π)

Another issue is when $x$$x$ is very large:

xxxxxxxxxx
ε = eps() # machine epsilon, 2^(-52)
x = 2*10.0^100
abs(sin(x) - sin(big(x)))  ≤  abs(sin(big(x))) * ε

But if we instead compute 10^100 using BigFloat we get a completely different answer that even has the wrong sign!

xxxxxxxxxx
x̃ = 2*big(10.0)^100
sin(x), sin(x̃)

This is because we commit an error on the order of roughly

when we round $2\ast {10}^{100}$$2*10^{100}$ to the nearest float.

Example (polynomial near root) For general functions we do not generally have relative accuracy. For example, consider a simple polynomial $1+4x+x^2$$1 + 4x + x^2$ which has a root at $\sqrt{3}-2$$\sqrt 3 - 2$. But

xxxxxxxxxx
f = x -> 1 + 4x + x^2
x = sqrt(3) - 2
abserr = abs(f(big(x)) - f(x))
relerr = abserr/abs(f(x))
abserr, relerr # very large relative error

We can see this in the error bound (note that $4x$$4x$ is exact for floating point numbers and adding $1$$1$ is exact for this particular $x$$x$):

Using a simple bound $|x|<1$$|x| < 1$ we get a (pessimistic) bound on the absolute error of $3{ϵ}_{\mathrm{m}}$$3 ϵ_{\rm m}$. Here f(x) itself is less than $2{ϵ}_{\mathrm{m}}$$2 ϵ_{\rm m}$ so this does not imply relative accuracy. (Of course, a bad upper bound is not the same as a proof of inaccuracy, but here we observe the inaccuracy in practice.)

5. High-precision floating-point numbers (advanced)

It is possible to set the precision of a floating-point number using the BigFloat type, which results from the usage of big when the result is not an integer. For example, here is an approximation of 1/3 accurate to 77 decimal digits:

xxxxxxxxxx
big(1)/3

Note we can set the rounding mode as in Float64, e.g., this gives (rigorous) bounds on 1/3:

xxxxxxxxxx
setrounding(BigFloat, RoundDown) do
  big(1)/3
end, setrounding(BigFloat, RoundUp) do
  big(1)/3
end

We can also increase the precision, e.g., this finds bounds on 1/3 accurate to more than 1000 decimal places:

xxxxxxxxxx
setprecision(4_000) do # 4000 bit precision
  setrounding(BigFloat, RoundDown) do
    big(1)/3
  end, setrounding(BigFloat, RoundUp) do
    big(1)/3
  end
end

In the problem sheet we shall see how this can be used to rigorously bound $\mathrm{e}$${\rm e}$, accurate to 1000 digits.

Differentiation

YouTube Lectures: Finite Differences Dual Numbers

We now get to our first computational problem: given a function, how can we approximate its derivative at a point? Before we begin, we must be clear what a "function" is. Consider three possible scenarios:

1. Black-box function: Consider a floating-point valued function $f^{\mathrm{F}\mathrm{P}}:D\to F$$f^{\rm FP} : D \rightarrow F$ where $D\subset F\equiv F_{\sigma ,Q,S}$$D \subset F \equiv F_{\sigma,Q,S}$ (e.g., we are given a double precision function that takes in a Float64 and returns another Float64) which we only know pointwise. This is the situation if we have a function that relies on a compiled C library, which composes floating point arithmetic operations. Since $F$$F$ is a discrete set such an $f^{\mathrm{F}\mathrm{P}}$$f^{\rm FP}$ cannot be differentiable in a rigorous way, therefore we need to assume that $f^{\mathrm{F}\mathrm{P}}$$f^{\rm FP}$ approximates a differentiable function $f$$f$ with controlled error in order to state anything precise.
2. Generic function: Consider a function that is a formula (or, equivalentally, a piece of code) that we can evaluate it on arbitrary types, including custom types that we create. An example is a polynomial: $p(x)=p_0+p_{1}x+\cdots +p_n{x}^n$$p(x) = p_0 + p_1 x + \cdots + p_n x^n$ which can be evaluated for $x$$x$ in the reals, complexes, or any other ring. More generally, if we have a function defined in Julia that does not call any C libraries it can be evaluated on different types. For analysis we typically consider both a differentiable function $f:D\to \mathbb{R}$$f : D \rightarrow {\mathbb R}$ for $D\subset \mathbb{R}$$D ⊂ {\mathbb R}$, which would be what one would have if we could evaluate a function exactly using real arithmetic, and $f^{\mathrm{F}\mathrm{P}}:D\cap F\to F$$f^{\rm FP} : D \cap F \rightarrow F$, which is what we actually compute when evaluating the function using floating point arithmetic.
3. Graph function: The function is built by composing different basic "kernels" with known differentiability properties. We won't consider this situation in this module, though it is the model used by Python machine learning toolbox's like PyTorch and TensorFlow.

We discuss the following techniques:

1. Finite-differences: Use the definition of a derivative that one learns in calculus to approximate its value. Unfortunately, the round-off errors of floating point arithmetic typically limit its accuracy.
2. Dual numbers (forward-mode automatic differentiation): Define a special type that when applied to a function computes its derivative. Mathematically, this uses dual numbers, which are analoguous to complex numbers.

Note there are other techniques for differentiation that we don't discuss:

1. Symbolic differentiation: A tree is built representing a formula which is differentiated using the product and chain rule.
2. Adjoints and back-propagation (reverse-mode automatic differentiation): This is similar to symbolic differentiation but automated, to build up a tape of operations that tracks interdependencies. It's outside the scope of this module but is computationally preferred for computing gradients of large dimensional functions which is critical in machine learning.
3. Interpolation and differentiation: We can also differentiate functions globally, that is, in an interval instead of only a single point, which will be discussed later in the module.

xxxxxxxxxx
using ColorBitstring

1. Finite-differences

The definition

tells us that

provided that $h$$h$ is sufficiently small.

It's important to note that approximation uses only the black-box notion of a function but to obtain bounds we need more.

If we know a bound on $f^{″}(x)$$f''(x)$ then Taylor's theorem tells us a precise bound:

Proposition The error in approximating the derivative using finite differences is

$$|f^{\prime }(x)-\frac{f(x+h)-f(x)}{h}|\le \frac{M}{2}h$$

where $M=\underset{x\le t\le x+h}{sup}|f^{″}(t)|$$M = \sup_{x \leq t \leq x+h} |f''(t)|$.

Proof Follows immediately from Taylor's theorem:

for some $x\le t\le x+h$$x ≤ t ≤ x+h$.

◼️

There are also alternative versions of finite differences. Leftside finite-differences:

and central differences:

Composing these approximations is useful for higher-order derivatives as we discuss in the problem sheet.

Note this is assuming real arithmetic, the answer is drastically different with floating point arithmetic.

Does finite-differences work with floating point arithmetic?

Let's try differentiating two simple polynomials $f(x)=1+x+x^2$$f(x) = 1 + x + x^2$ and $g(x)=1+x/3+x^2$$g(x) = 1 + x/3 + x^2$ by applying the finite-difference approximation to their floating point implementations $f^{\mathrm{F}\mathrm{P}}$$f^{\rm FP}$ and $g^{\mathrm{F}\mathrm{P}}$$g^{\rm FP}$:

xxxxxxxxxx
f = x -> 1 + x + x^2     # we treat f and g as black-boxs
g = x -> 1 + x/3 + x^2
h = 0.000001
(f(h)-f(0))/h, (g(h)-g(0))/h

Both seem to roughly approximate the true derivatives ($1$$1$ and $1/3$$1/3$). We can do a plot to see how fast the error goes down as we let $h$$h$ become small.

xxxxxxxxxx
using Plots
h = 2.0 .^ (0:-1:-60)  # [1,1/2,1/4,…]
nanabs = x -> iszero(x) ? NaN : abs(x) # avoid 0's in log scale plot
plot(nanabs.((f.(h) .- f(0)) ./ h .- 1); yscale=:log10, title="convergence of derivatives, h = 2^(-n)", label="f", legend=:bottomleft)
plot!(abs.((g.(h) .- g(0)) ./ h .- 1/3); yscale=:log10, label = "g")

In the case of $f$$f$ it is a success: we approximate the true derivative exactly provided we take $h=2^{-n}$$h = 2^{-n}$ for $26<n\le 52$$26 < n \leq 52$. But for $g$$g$ it is a huge failure: the approximation starts to converge, but then diverges exponentially fast, before levelling off!

It is clear that $f$$f$ is extremely special. Most functions will behave like $g$$g$, and had we not taken $h$$h$ to be a power of two we also see divergence for differentiating $f$$f$:

xxxxxxxxxx
h = 10.0 .^ (0:-1:-16)  # [1,1/10,1/100,…]
plot(abs.((f.(h) .- f(0)) ./ h .- 1); yscale=:log10, title="convergence of derivatives, h = 10^(-n)", label="f", legend=:bottomleft)
plot!(abs.((g.(h) .- g(0)) ./ h .- 1/3); yscale=:log10, label = "g")

For these two simple examples, we can understand why we see very different behaviour.

Example (convergence(?) of finite difference) Consider differentiating $f(x)=1+x+x^2$$f(x) = 1 + x + x^2$ at 0 with $h=2^{-n}$$h = 2^{-n}$. We consider 3 different cases with different behaviour, where $S$$S$ is the number of significand bits:

1. $0\le n\le S/2$$0 ≤ n ≤ S/2$
2. $S/2<n\le S$$S/2 < n ≤ S$
3. $S\le n$$S ≤ n$

Note that $f^{\mathrm{F}\mathrm{P}}(0)=f(0)=1$$f^{\rm FP}(0) = f(0) = 1$. Thus we wish to understand the error in approximating $f^{\prime }(0)=1$$f'(0) = 1$ by

Case 1 ($0\le n\le S/2$$0 ≤ n ≤ S/2$): note that $f^{\mathrm{F}\mathrm{P}}(h)=f(h)=1+2^{-n}+2^{-2n}$$f^{\rm FP}(h) = f(h) = 1 + 2^{-n} + 2^{-2n}$ as each computation is precisely a floating point number (hence no rounding). We can see this in half-precision, with $n=3$$n = 3$ we have a 1 in the 3rd and 6th decimal place:

xxxxxxxxxx
S = 10 # 10 significant bits
n = 3  # 3 ≤ S/2 = 5
h = Float16(2)^(-n)
printbits(f(h))

Subtracting 1 and dividing by $h$$h$ will also be exact, hence we get

which shows exponential convergence.

Case 2 ($S/2<n\le S$$S/2 < n ≤ S$): Now we have (using round-to-nearest)

Then

We have actually performed better than true real arithmetic and converged without a limit!

Case 3 ($n>S$$n > S$): If we take $n$$n$ too large, then $1\oplus h=1$$1 ⊕ h = 1$ and we have $f^{\mathrm{F}\mathrm{P}}(h)=1$$f^{\rm FP}(h) = 1$, that is and

Example (divergence of finite difference) Consider differentiating $g(x)=1+x/3+x^2$$g(x) = 1 + x/3 + x^2$ at 0 with $h=2^{-n}$$h = 2^{-n}$ and assume $n$$n$ is even for simplicity and consider half-precision with $S=10$$S = 10$. Note that $g^{\mathrm{F}\mathrm{P}}(0)=g(0)=1$$g^{\rm FP}(0) = g(0) = 1$. Recall